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Abstract 

The extended Finite Element Method (XFEM) is used to solve interface problems with 
an unfitted mesh. We present an implementation of the XFEM in the FEM-library deal.IF 
The main parts of the implementation are (i) the appropriate quadrature rule; (ii) the 
shape functions for the extended part of the finite element formulation; (iii) the boundary 
and interface conditions. We show how to handle the XFEM formulation providing a code 
that demonstrates the solution of two exemplary interface problems for a strong and a 
weak discontinuity respectively. In the weak discontinuity case, the loss of conformity due 
to the blending effect and its remedy are discussed. Furthermore, the optimal convergence 
of the presented unfitted method is numerically verified. 


1 Introduction 

The extended Finite Element Method (XFEM) is a flexible numerical approach developed for 
general interface problems. Numerical methods to solve interface problems can be classihed 
as fitted or unfitted methods. In the hrst case, the methods use a fitted mesh approach such that 
the interface is composed of element sides. The generation of a htted mesh in case of complex 
interface geometry can be very time consuming. In many cases it cannot be done without 
handwork using a program for mesh generation. In the unhtted case, the mesh is independent 
of the interface position and therefore unhtted methods are highly hexible. Since standard 
hnite element methods perform poorly in the unhtted case, diherent alternative approaches 
have been introduced in the last years. 

The XFEM is a partition of unity hnite element method (PUFEM). The hrst formulation 
of the PUFEM has been derived in the work of Melenk and Babuska [13]. The main features 
of this method can be summarized as follows: 

• a priori knowledge about the local behavior of the solution can be included in the formu¬ 
lation; 

• arbitrary regularity of the FE spaces can be constructed; 

• the approach can be understood as a meshless method; 

• it is a generalization of the h, p and hp version of the FEM. 
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In particular two important aspects are essentially relevant for this method: local approxima- 
bility and the capability to enforce inter-element continuity, i.e. conformity. Among different 
PUFEM approaches, the generalized hnite element method (GFEM) and the extended hnite 
element method are the most versatile and the most used in many applications. Their elabo¬ 
ration developed from the area of meshfree methods [7] and are based on the same principles: 
partition of unit and degree of freedom enrichment [10]. An overview of these methods can be 
found for example in [1, 2]. 

The XFEM strategy to solve a problem with strong or weak discontinuities, i. e. disconti¬ 
nuities of the solution or of the fluxes respectively, is to extend the approximation space with 
discontinuous basis functions or basis functions with a kink, respectively. Since the discontinu¬ 
ities are typically local features, the XFEM offers great flexibility by using a local modihcation 
of the standard FEM methods. In fact, it avoids the use of complex meshing, which is substi¬ 
tuted by a specihc distribution of the additional degrees of freedom (DoFs). 

The XFEM has broad use in different disciplines including fracture mechanics, large defor¬ 
mation, plasticity, multiphase flow, hydraulic fracturing and contact problems [12]. However, 
the hrst developments of the XFEM were done to simulate crack propagation [15]. Further 
applications for the XFEM in material science comprise: problems with complex geometries, 
evolution of dislocations, modeling of grain boundaries, evolution of phase boundaries, model¬ 
ing of inclusions and homogenization problems. In particular, the combination of the XFEM 
with a level set approach [18] has been shown to be a very versatile tool to solve the above 
class of problems. In the level set approach two strategies can be used to dehne the interfaces: 
(i) an analytical description of interfaces as the iso-zero of a function can be given or (ii) data 
from an image segmentation can be used to dehne interfaces locally or globally. 

The goal of this work is to present all essential steps for the implementation of the XFEM. In 
addition, we make available a code (contact the authors to get a copy of it) that can be further 
extended for specihc applications. The practical implementation is done in the open source FEM 
library deal.II [3]. As application we consider an interface problem with a weak or a strong 
discontinuity. In the case of a strong discontinuity we consider only weak interface conditions of 
Robin type. In particular, we do not consider Dirichlet conditions on the interface. In this case 
the formulation of the problem has to be changed and possible formulations are based either 
on the Nitsche’s method [17], see for example [11, 5], or on a Lagrange multipliers method 
[14, 4]. The extension of the code including the Dirichlet case is left for a further development. 
In addition, we illustrate the problem of the blending effect, i.e. of the loss of conformity in 
the elements adjacent to the interface, and show the implementation of a standard method to 
restore the full convergence behavior. The focus of this work is on implementation aspects for 
the stationary XFEM. We do not consider therefore moving interfaces. The development of 
a proper time stepping technique and an adequate quadrature rule goes beyond the scope of 
this paper. We present examples in the two-dimensional case. Specihc aspects related to the 
extension to three-dimensional problems are also discussed. 

This article is organized as follows. In section 2 we introduce the general interface problem 
and the level set method. In section 3 the formulation of the XFEM for strong and weak 
discontinuities is depicted. Furthermore, the blending effect in case of weak discontinuities is 
discussed. In section 4, we briehy report some theoretical results on the convergence of the 
XFEM. In section 5 we describe in detail our implementation in deal.II. Specihcally we discuss 
the XFEM quadrature rule and the application of boundary conditions to cut cells. In the hnal 
section 6 we show some numerical results on problems with weak and strong discontinuities 
including convergence tests. 
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2 Interface problems 

2.1 Problem setting 

Let be a bounded domain in with boundary dQ. The considered model problem is the 
stationary heat conduction. We construct an interface problem by taking a domain hi divided 
in two sub-domains hli and ^2 by a line T, called interface. 

We consider three cases of interface problem with discontinuity on T. Case I: the solution 
has a discontinuity Qs] Case II: the solution has a weak discontinuity Case III: the solution 
has a discontinuity which depends on some functions gi, g 2 as shown below. The problem can 
be formulated as 

Problem 2.1 (Interface problem). Given the function f, the strong and weak discontinuity at 
T, gs and g^ respectively, or the functions gi, g 2 , find the solution u of the following system 



-V ■ (/iiVui) 

= f 

inDj, 

(la) 


Ui 

= 9, 

on dfli n D, 

(lb) 

Gase I : 

[u] 

= 9 s 

onT, 

(Ic) 

Gase II : 

[h'Vu ■ n] 

9w 

onT, 

(Id) 


[u] 

= 0 

onT, 

(le) 

Gase III - 

• Ui 

= 9i{'ai,U2) 

onT, 

(If) 


fori = 1,2, with [u] =ui — U 2 and [fvVu-n] = fii'Vui-ni + g, 2 '^U 2 -n 2 , where Ui is the restriction 
of u to fli, Hi is a constant assumed positive, and nt is the outward pointing normal to 12^ at T, 
see Figure 1. 

We consider a general discretization with finite elements and use the notation with subscript 
h to indicate discretized functions. 


Problem 2.2 (Discrete interface problem). Given the same data as the continuous problem 
above, find the solution Uh of the following discretized system 


-V ■ (/XiVuft,*) = /, 


Gase I: [uh] = gs 

Gase II : [h'^Uh ■ n] = g^, 

[uh] = 0 

Gase III. j • Uj 


inDj, 

(2a) 

oncIDj n D, 

(2b) 

onT, 

(2c) 

onT, 

(2d) 

onT, 

(2e) 

onT, 

(2f) 


for i = 1,2, with [uh] = Uh,i - Uh ,2 and ■ n] = ■ ni + H 2 ^Uh ,2 ■ n 2 , where Uh,i is 

the restriction of Uh to D,. 

Note that the same interface T of the continuous problem is used also for the discretized 
problem if an exact quadrature formula can be employed. An exact representation of T used 
in the discretized problem can be obtained adopting the level set method [18]. 
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2.2 Level set method 


The level set method is used to implicitly define the position of the interface T independently 
of the underlying mesh used to discretize the interface problem. The interface is defined by 
a scalar function 0 : hi —)■ M that is (uniquely) zero on T and has different sign on different 
sub-domains: 


0 = 0 on T, 

0 < 0 in hli, (3) 

0 > 0 in r22- 

Typically, the distance function (with sign) to the 

(bix) = ± min I 

yer ' 

This is not the only choice, but it is often used because it can be exploited to calculate the 
normal vector at any point on T, since V0/| V0| represents the normal vector if 0 is the distance 
function to T. Following a standard definition of XFEM we use the level set function to extend 


interface F is used as level set function 
\x-y\\. 



the space of finite elements with functions that incorporate the needed discontinuity. 

3 Extended finite elements 

Generally, Galerkin finite elements are defined as the triple 


(r,gp,s), (4) 

where T is the mesh, Qp is the space of test and trial functions and S is a set of linear functionals 
that defines the degrees of freedom of the FEM formulation. To build the XFEM in deal.II, 
we consider an extension of Lagrange finite elements. These are FE for which the degrees of 
freedom are the values of the test functions at the nodes of the mesh elements. Since our 
implementation is done in deal.II we consider only quadrilateral elements iL G T. 

The space Qp is the space of polynomial functions of degree at most p 

Qp = {fix) = 

ai,...,ad<p 
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defined on a nnit cell K = (0, l)*^. The test and trial fnnctions on a real cell K are obtained 
throngh a transformation a : K ^ K of a. fnnction from Qp. We will use the notation G Qp 
to indicate that the transformation of the function onto the unit cell belongs to Qp, i.e., 
(j~^{(Pk) e Qp. 

The scope of this work is to describe the implementation of the XFEM in deal.II, so we 
restrict our test cases to linear problems without loss of generality. We consider a general 
elliptic bilinear form and a linear functional: 

a : V xV (5) 

( 6 ) 

where V is an appropriate Hilbert space. The general weak formulation of our test cases is: 

Problem 3.1. Find u ^ V such that 

a{u,F) = f{F), V(peH. (7) 

A typical choice for V is the Hilbert space H^{Q) where H is the problem domain. 

The discrete approximation of problem 3.1 using hnite elements is 

Problem 3.2. Find Uh G 14 such that 

a{uh,(ph) = fiFh), V(p/i G 14, (8) 

where 14 is the hnite dimensional space of if^-conform functions 

W = {(P G 1/ : (p4 G Qp WK G T}, 

The solution vector is a linear combination of the basis functions of 14 

n 

Mx) = '^UjNj{x). (9) 

j=i 

In the following we restrict our formulation to the space of bilinear functions Qi, i.e. the shape 
functions W are piecewise bilinear and globally continuous. 

For interface problems of the type 2.1 there is the need to approximate a solution with 
a discontinuity. As illustrated above in the interface problem 2.1, we consider three cases of 
boundary conditions on the interface leading to a weak discontinuity or a strong discontinuity. 
In case of a weak discontinuity the standard Qi space can reach the best convergence rate (for 
a solution smooth enough) only if the mesh is htted with the weak discontinuity. In case of 
strong discontinuity, we consider convergence in a norm in the space F[^{Qi U 122), since in the 
space 77^(12), where 12 = 12i U 122 U F, the solution is not continuous and therefore it does not 
belong to 77^. Also in this case, the standard Qi space can achieve the best convergence rate 
(in 77^(12i U I 22 )) only if the degrees of freedom he on F. 

If the interface cuts some elements K, a better approximation can be achieved by incorpo¬ 
rating the discontinuity in the space in which we approximate the solution. 


In the considered XFEM formulation we extend therefore the space Qi with some additional 
shape functions that represent the given discontinuity. This is obtained by enriching the degrees 
of freedom of the elements cut by the interface. We will use the notation F for the standard 
degrees of freedom and I* for the set of the extended degrees of freedom, respectively represented 
with single points and double points in Figure 2. The set of all degrees of freedom is denoted 
7. 

In the next subsection we construct the XFEM shape functions. 
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Figure 2: Single dots depicts normal degrees of freedom. Double dots depicts the extended 
degrees of freedom. 


3.1 Strong and weak discontinuity 

In case of strong discontinuity, a typical function with a jump along F is the sign function: 

sign : M —)■ {—1,0,1} 

{ 1 for a; > 0 

0 for a; = 0 

— 1 for a; < 0. 

Since the jump is at the point a; = 0, the sign function applied to the level set function can be 
used to obtain a function with a jump along the interface. 

In case of weak discontinuity, a function with a kink can be used, as for example the absolute 
value: 


abs : 


-)■ 


M+ 


I X 

for 

0 

for 

1 —X 

for 


a; > 0 
a; = 0 
a: < 0, 


which applied to the level set function dehnes a weak discontinuity along the interface. 
In the following we use the general notation 


■^ : D —>■ 
%Ij{x) = 


fsign(0(a;)) 
I abs(0(a:)) 


for strong discontinuity 
for weak discontinuity. 


is called enrichment function. 

In Figure 2 the extended degrees of freedom are depicted. These are additional Lagrangian 
degrees of freedom dehned on a subset of existing mesh nodes. To construct the XFEM, 
additional shape functions have to be dehned. Following the above construction, we take 
functions that have a discontinuity along F 


Mi{x) := Ni{x)tlj{x). 

We consider thus the following discrete spaces for test and trial functions 


( 10 ) 


6 























strong discontinuity 


: (p\k e Qi, ^\k[ e Qi,i — 1,2}, 

where K are standard cells and K' are the cells cut by the interface and K[ are the two 
parts of the cut cell K' that have the interface as common edge. 

• weak discontinuity: 


Vh '■ ^\k G Qi) v\k' G Qi © |0|Qi}. 


The discrete solution is now defined using the enriched basis 

Uh{x) = '^UiNi{x) + ^ ajMj{x). 
ier jei* 


From the practical point of view, it is desired that the discrete solution has the Kronecker delta 
property 

Uh{xi) =Ui i = 1,... ,n, (11) 

where Xi is the position of the degree of freedom. To obtain this property, the extended 
basis functions are shifted, i.e. we use a different basis. The modihed extended basis function 
of the node i becomes 

Mi{x) = Ni{x){'4){x) - '4){xi)). ( 12 ) 

In Figure 3 two extended basis functions are depicted, for the cases of weak and strong discon¬ 
tinuity. The extended basis functions depend on the position of the interface. In case of weak 




Figure 3: XFEM basis functions for weak discontinuity (left) and strong discontinuity (right). 


discontinuity, the use of standard enriched basis functions can lead to the loss of conformity in 
the elements adjacent to the cut cells. This so called blending effect introduces a reduction of 
the convergence rate as it is shown later in the numerical experiments. 


3.2 Blending effect 

In this subsection we discuss the blending effect. Let’s consider the basis functions for the weak 
discontinuity 


Mi{x) = A^i(a;)(abs(())(a;)) - abs(())(xi))). 


(13) 
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which is given by the product of a polynomial basis function Ni and the level set function, 
which is in general a function. Considering the term abs(0(a:)) — ahs{(j){xi)) along an edge 
E of the cell where the function is not zero, it is 


abs(0(a;)) — abs(0(a;j)) 


= constant 
7 ^ constant 


for E II r 
for E\{V. 


(14) 


Since the function Ni along the edge E is linear, the product with the enrichment function 
gives a function that is non linear. In Figure 4 an extended basis function is depicted that 
shows the blending effect. The depicted function has a nonlinear behavior on the right side 



along the x axis and on the left side along the y axis. The nonlinearity along the x axis is 
not a problem, since the edge is in common with a cut neighbor cell, which is also enriched 
in the same manner. On the other side, the edge along the y axis causes problems, because 
the neighbor element is not cut by the interface and has only standard basis functions, which 
are linear on the common edge. We have thus linear behavior on one side and nonlinear on 
the other. Due to the discontinuity along this edge the enriched space Vj^ is not if^-conform 
anymore. Only in the case that the interface is parallel to the edges there is no blending effect, 
because the enrichment function is constant along such edges, see (14). 

There are different approaches to recover the if^-conformity: 

• use of higher order elements in the standard space 14 [8] , 

• smoothing techniques as used in [19], 

• use of a corrected XFEM formulation adding a ramp function [9]. 

Note that the use of higher order elements works only if the extended functions on the cell 
edges have polynomial behavior. In this work we use the third method using a correction with 
the ramp function: 


’’(a^) := (15) 

iei° 

where 1° depicts the set of standard degrees of freedom that lie on cut cells. 

The idea is to enrich not only the cells that are cut by the interface, but also the neighbor 
cells that are called blending cells. In these cells the basis functions are modified so that they 
are nonlinear on the common edge with a cut cell and linear on the other edges recovering the 
global continuity. 
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Figure 5: Ramp function in a blending cell which neighbor cell on the right is a cut cell. 


Figure 5 shows the ramp function on the neighbor of the cut cell depicted in Figure 4. The 
ramp function has the value 1 along the edge that creates the blending effect and is zero on the 
opposite edge. By multiplying the XFEM basis functions with the ramp function, new basis 
functions are defined that impose the continuity along the common edge with cut cells thus 
deleting the blending effect. Therefore, in the blending cells we use the extended functions: 

Mi{x) := Ni{x){aJos{4>{x)) — ahs{4>{xi)))r{x). (16) 

which are depicted in Figure 6. It can be observed that these additional functions behave 



Figure 6: Extended basis functions on blending cells. 


nonlinearly along the common edge with a cut cell and a blending cell (in Figure 6 depicted 
on the right side and on the bottom side respectively) and go to zero on common edges with 
normal cells (left and top side). 

4 A note on a priori error estimation 

In this section we briefly present some known results on a priori convergence estimates for the 
XFEM for interface problems. A result on optimal convergence rate of the XFEM for crack 
propagation can be found in [16]. Let’s consider the interface problem of the type III using the 
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same notation as in Problem 2.1 


-piAu = f 

in 12j 

(17a) 

u = 0 

on 512 

(17b) 

PldniUi = Q-iiUi + 0:12^2 

on F 

(17c) 

h'25n2'W2 = Oi2lUi + CX 22 U 2 

on F. 

(17d) 


Note that continuity along the interface P is not enforced. For this problem we consider the 
weak formulation 

Problem 4.1 (Interface problem: Weak formulation). Let the problem data be regular enough 
and let hlj be two domains with smooth boundaries so that the regularity u G U 122) is 

assured. Find u E V, such that for all (p E V it is 

a{u, (p) = n,un2 (18) 

with 

a{u, (p) = (pi Vm, Vv 7 )oi + (/i2Vu, Vp)n2 

— (OiiMi, (^i)r — (012^2, 9 ^l)r — («21 Mi, P2)v — («22'U2, </ 52 )r, 

and V = hfg (12i U 122; 512), where we have used the notation 512 := 5(12i U I 22 ) \ P. 

Let’s consider the mesh T and the finite dimensional space Vh <ZV: 

Vh := {p>h G V : (phlxi G Qi, Fh\Ki = 0 V 'Ph\K2 = 0) eT, i = 1,2} (20) 

with Ki := il n l2j. The basis functions ph are the unhtted basis functions used in [11] 
to show the convergence results. It can be shown that XFEM basis functions together with 
standard basis functions build a basis for the hnite dimensional space I 4 . In fact, as observed 
by Belytschko in [6] the unhtted basis functions in [11] are equivalent to the XFEM basis 
functions. Therefore it is 


Vh = spanjW, Mj : i E j E /*}, (21) 

where /' and I* are dehned as in section 3. We consider the following XFEM approximation 
of the above problem 

Problem 4.2 (Interface problem: Discrete formulation). With the same data as the above 
continuous problem, find Uh E Vh, such that for all ph £ 14 it is 

a{uh,Ph) = if:Fh)n, (22) 


with 12 = 12i U 122 U F. 

Under these conditions following [11] it can be shown that for the hnite element solution Uh 
using the XFEM it is 


and 


1|V('U — 'Uh)||oiu02 < ch||-u||//2(QjUQ2)) 


(23) 


\u 


m / i || oiu 02 4 chf\\u\\m(n^un2)- 


(24) 
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Note that to show these error estimations using the results in [11], one has to perform a change 
of basis since Hansbo and Hansbo use different basis functions than the XFEM ones as pointed 
out above. 

Therefore, full convergence behavior has to be expected in our numerical tests with strong 
discontinuity. Indeed as it will be shown later, also in case of weak discontinuity, using the 
ramp correction on blending cells, we observe the same convergence behavior. Nevertheless, the 
a priori convergence estimates for the weak discontinuity case cannot be derived in the same 
way following the work of Hansbo and Hansbo. 


5 Implementation in deal.II 

To implement the XFEM in deal.H we consider the scalar interface problem 2.1 as a vectorial 
problem. The solution is therefore represented with two components. One component is the 
standard part of the solution and the other is the XFEM extension. The standard part of 
the solution exists in all cells, whereas the extended part exists only in the cells that are cut 
by the interface, and their neighbors (blending cells) in case of weak discontinuity. Due to 
an implementation constraint, both components of the vector-valued function must however 
exist in all cells. Therefore, since the degrees of freedom (in particular those belonging to the 
extended part) are distributed to all cells, we have to extend with the zero function the part of 
the extended solution over uncut cells. By solving the system of equations it must be ensured 
that the zero extension of the solution remains also identically zero. 

This is obtained in the implementation in deal.H with two main objects: the class FeNothing 
and the class hp::DoFHandler. The object FeNothing is a hnite element class that has zero 
degrees of freedom. The object hp::DoFHandler allows to distribute different hnite element 
types on different cells. Since we use the vector-valued hnite element FESystem with two com¬ 
ponents, we can arbitrarily assign to each cell the two types of hnite element either FeNothing 
or Qi. We assign thus a Qi hnite element object to the hrst component (standard FE) of all 
cells, while we consider the following three cases for the second component of the FESystem: 

(i) for cells cut by the interface we use a Qi object to dehne the extended part of the space; 

(ii) for the blending cells we use a Qi object to dehne the extended part of the space including 
the ramp functions. 

(iii) for rest of the cells we use a FeNothing object (the extended part is set to zero); 

The distribution of degrees of freedom with hp::DoFHandler is controlled by the value 
active-fe-index of the cell iterator. 

5.1 Quadrature formula 

An essential part of the XFEM implementation is the quadrature formula. Generally, a quadra¬ 
ture formula such as the Gauss quadrature designed to integrate smooth functions would fail to 
integrate a function with a discontinuity. A proper quadrature formula must take into account 
the position of the interface to integrate a function with a weak or strong discontinuity. In our 
deal.H implementation this is done by subdividing the cut cells in sub-elements, on which a 
standard quadrature formula can be used. Since we restrict our implementation to standard 
elements in deal.H, i.e. quadrilateral elements in 2D, the subdivision is done by quadrilateral 
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sub-elements. In 2D there are only four types of subdivisions and the respective rotated vari¬ 
ants, see Figure 7, i.e. four of type (a), eight of type (b), two of type (c) and two of type 
(d). Let’s use the notation K for the cell in real coordinates, K and K for the unit cell, and 





Figure 7: Subdivisions of the unit cell. 

Si,... ,Sn for the subcells of K. Since in deal.II the quadrature formula is dehned for the refer¬ 
ence unit cell, we have to construct a quadrature formula with points and weights to integrate 
a transformed function on the unit cell. To this aim, the cut cell K is transformed into the unit 
“cut” cell K and with the help of the transformed level set function it is subdivided according to 
the schemes in Figure 7. Each subcell is then transformed into the unit “uncut” cell K in order 
to calculate the local quadrature points and weights through the standard tools in deal.II. Note 
that to make clear the use of the unit cell in different situations we use two notations for it, i.e. 
K and K. Subsequently they are transformed back to the cell in “real” coordinates that in this 
case are the coordinates of the unit cell K. Different transformations are used to transform K 
to K and each Si into K as depicted in the same hgure. The transformation a: K ^ K is 
the standard transformation used in deal.II. Furthermore, to construct the appropriate XFEM 
quadrature formula we transform the subcell, see Figure 8, through the transformations: 

ai,..., an ■ K ^ Si ,..., Sn- (25) 

The XFEM quadrature formula is therefore derived by a standard quadrature formula. 


K k k 



Figure 8: Transformations to construct the XFEM quadrature formula. 


For given quadrature points Xi and weights Wi of a standard formula, i = 1,... ,m, the 
XFEM quadrature on k is dehned through the points 

yij = akxi), i = l,...,m, j = l,...,n, (26) 
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and weights 


Wij = Wi det{Vaj{xi}). (27) 

We define the degree of exactness of a qnadratnre formnla as the maximal degree of polynomial 
fnnctions that can be exactly integrated on an arbitrary domain. The optimal position of the 
qnadratnre points depends on the shape of the domain on which the integration is done. The 
degree of exactness of the XFEM qnadratnre formnla is at least of the same degree as the 
standard formnla from which it is derived. In addition, it allows by constrnction to integrate 
discontinnous fnnctions along the interface. From the implementation point of view this formula 
is highly flexible because it is built as a combination of standard formulas. The construction 
of an XFEM quadrature formula is simplified in our implementation since the subdivision of a 
cut cell is done by using the same type of cells, i.e. quadrilateral cells in our two-dimensional 
case. 

We would like to underline, that the XFEM formula is not optimal from the theoretical point 
of view, since it does not use the minimal number of points needed to integrate a given function 
over the specific subdivisions. In fact, the degree of exactness of the formula could be higher 
than the one inherited by the standard formula, but never less accurate. Therefore, the use 
of the XFEM formula can result in unnecessary higher costs to integrate a given function and 
therefore higher costs, for example, in assembling system matrices and vectors. As an example, 
let’s consider the integration of a quadratic function on the triangle resulting from the cut 
depicted in Figure 7 in the cases (a), (b) and (d). The unit cut cell is divided in two parts, one 
triangle and one pentagon. For the pentagon part, quadrature rules would be necessary that are 
not typical on finite elements codes and therefore are not present as standard implementation. 
Therefore the division in two parts of the pentagon to obtain two quadrilaterals on which we 
can use standard formula is a desired feature. On the contrary on the triangle part one could 
use a more efficient formula, for example we could use a symmetric Gauss quadrature with 3 
points, while the XFEM quadrature rule is build with 12 points as depicted in Figure 9. The 
code can be slightly more efficient changing the quadrature rule for the triangles obtained by 
the cut. We do not consider this modification for two reasons. On one side, we expect a great 
gain in computing time only in cases with a very high number of cut cells. In a two dimensional 
numerical test not shown here we have observed about 13% reduction of computing time using 
a three point quadrature formula for the triangle subdivisions instead of the XFEM one. On 
the other, we want to produce a code that can be used in a dimension independent way. 




(a) (b) 

Figure 9: Sketch of the position in a triangle of the XFEM quadrature points (a) and of the 
symmetric Gauss formula (b). 

In the two dimensional case the two parts of a cut cell can only be quadrilaterals, triangles 
or pentagons. In the three dimensional case there are 15 cases (and their respective symmetric 
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configurations) and the partitions are more complex polyhedra as can be seen in the two 
depicted cases of Figure 10. Therefore, in our implementation instead of considering a complex 
quadrature formula that can cope with all possible subdivisions, we apply subdivisions using 
only standard cells of deal.II in two and three dimensions, i.e. quadrilaterals and hexahedra 
respectively. The extension to the three dimensional case is theoretically straightforward. In 
practice the subdivision of hexahedral cut cells in hexahedral subcell is not a trivial task and 
it is left for a forthcoming work in which we will consider a comparison of different quadrature 
rules. 



5.2 Boundary and interface conditions 

This section is dedicated to the boundary and interface conditions. We consider Dirichlet 
and Neumann boundary conditions on the external boundary. Furthermore, we describe the 
implementation of Robin interface conditions on the interface F. 


5.3 Dirichlet and Neumann boundary conditions 

In the following we consider the Dirichlet boundary condition (lb) on the boundary of the 
domain. Nevertheless, the case with Neumann conditions can be treated in a similar way. 

The Dirichlet boundary condition in deal.II is set by an appropriate modihcation of the sys¬ 
tem matrix and right hand side, and optionally performing one step of the Gaussian elimination 
process to recover original properties of the matrix as, e. g., the symmetry. 

Owing to the Kronecker-delta property of the XFEM formulation, no special care has to be 
taken in case the Dirichlet condition is given by a continuous function. The degrees of freedom 
associated with the standard part of the FE are used to set the boundary values, while the 
degrees of freedom of the extension are set to zero at the boundary. In case of discontinuous 
boundary condition, see Figure 11 on the left side, the extended FE are used to approximate 
the discontinuity. On the edge at the boundary there are two extended degrees of freedom, 
whose basis function are discontinuous at one point. Due to the shift (13) used to construct 
the extended basis functions, these are nonzero only on one side of the edge. Therefore they 
can be used uncoupled to approximate the discontinuous Dirichlet value. 

Let’s consider the point Xc of the discontinuity of g on the edge with vertex xi and 0 : 2 , and 
the two limit values of the Dirichlet function gi{xc) and g 2 {xc). The two degrees of freedom oi 
and 02 of the extended part at the boundary are uniquely dehned by the following limit; 

lim Uh{x) = gi{xc), (28) 

X—^'Xc, XGofli 
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Figure 11: Discontinuous boundary condition. 


which leads to 


Oi = 


02 = 


g2{xc) - gi{xi)Ni{xc) - g2{x2)N2{xc) 
2Ni{xc) 

gijxc) - gi{xi)Ni{xc) - g2{x2)N2{xc) 
-2N2{Xc) 


(29) 

(30) 


In case of a weak discontinuity at the point Xc, the formulation needs a similar condition as 
in (28) for the normal derivatives. The two extended degrees of freedom are therefore coupled 
and their value can be calculated solving a system of dimension 2x2. 


5.4 Robin interface conditions 


Following the notation of (2f) (strong discontinuity), to simplify the description of the interface 
conditions we assume that gi and g 2 are linear in both arguments. For the variational formula¬ 
tion of the discrete interface problem 2.2, with interface conditions (2f), we dehne the bilinear 
form 

^h) = a{iph, ^IJh)n + a{(ph, '0h)r, 


with 


a{iph,'iph)n = (/iiV(p/i, V^/'h)ni + {g^2'^<Ph,'^'iph)n2 


and the boundary integrals 


a{iph,'il^h)r = {gi{,g^h,i,g^h, 2 ),'il^h,i)r + (fi'2(v5h,i,</3fe,2),^h,2)r, 

where the subscript 1 and 2 denotes the limit to the interface of the restriction on the subcells 
of the basis and test functions, e.g. for a; G F 


= lim iph{x) (31) 

The bilinear form is used to build the system matrix and the scalar product has to be im¬ 
plemented considering all mixed products between standard and extended part of the hnite 
element formulation. 

In the case of continuous basis functions (standard FE part) y:>h,i and (fh ,2 coincide. On 
the contrary, in the extended part of the XFEM formulation the basis and test functions are 
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Figure 12: Exact solution of problem 6.2. 


discontinuous and one of the two functions in the scalar product vanishes due to the shift (12). 
To determine the value of the limit (31) we use the systerri-to-component Jndex in deal.II, which 
is a pair containing the component of the current DoF and the index of the shape function of 
the current DoF. With the help of the level set function we can determine on which side of the 
interface the current DoF lies. This uniquely determines the zero part of the function. 


6 Numerical examples 

In this section we consider the solution of a linear elliptic interface problem in case of strong 
and weak discontinuity. Furthermore, we show the effect of the blending cells To simplify the 
notation in the following we use the symbol D instead of VL^. The latter is the approximation 
of the boundary given by the adopted hnite element formulation. 

6.1 Weak discontinuity 

Let’s consider the following domains Di := {x G : ||x ||2 < 0.5}, D 2 := {x G : 0.5 < 

||x ||2 < 1} and D := Di U 122. with the interface F = {x G : ||x ||2 = 0.5}. We dehne the 
following problem 

Problem 6.1 (Weak discontinuity). Given pi = 20 and fi 2 = 1, find the solution u 


-V ■ (/UjVu) = 1 

in Dj, 

(32) 

u = 0 

on 912, 

(33) 

H = o 

on F, 

(34) 

[lldnU] ^ 0 

on F. 

(35) 

The exact solution of problem 6.1 is the function 





(36) 

( J- .(-l . WxP + D.'i 

U-a-Nn 

X G 

X G 122. 



In Figure 12 a XFEM approximation of this problem is shown. 
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Problem 6.2 (Weak formulation of weak discontinuity). Find u G Hq{Q), so that 

if^Vu,V^)n = {f,Fh (37) 

with ^ = III in fli and fi = fi 2 in 0.2- The conditions (33) - (34) are naturally fulfilled by the 
weak formulation. 

We consider in the following two approximations calculated using two different finite di¬ 
mensional spaces to show the blending effect. For a given mesh T let’s consider the three 
subsets 

• Tout '■= {K gT iiCnr 7 ^ 0 }: the set of cells cut by the interface; 

• Tbi- the set of blending cells, i.e. those cells that are neighbors of cut cells; 

• Tstd ■= {K ^ T K (fi (Tbi U Tstd)}'- the set of standard cells. 

We can thereby dehne the ifi-conform space 14 and the non-conform space 14: 

14 := {Fh -FiiIk e Qi for K G (Tstd U Ti), 

Fh\K, G Qi © \(f)\Qi for K G Tut, 
cpbeC(n\T)} 

14 := {Fh e Hq(Q) -.iphlK e Qi for K G Tstd, 

Fhlx, G Qi © |0|Qi for K G Tut, 

Fhlx e Qi ®r\(j)\Qi for K G Ti, 

^Ph G ^(f])}. 

With this notation the discretized problem is given by: 

Problem 6.3 (Discrete formulation of weak discontinuity). With the data from Problem 6.2, 
find Uh G 14; so that 

(fiTuh,Tiph)n = (f,Fh)n M^Ph^Vh, (38) 

with either Vh = Vh or Vh = Vh. 

We use an unhtted mesh, i.e. the interface F intersects some cells. Therefore the underlying 
computing mesh is a shape regular mesh as shown in Figure 13 (a). Figures 13 (b) and (c) show 
the subdivisions in subcells for two level of the mesh. Remind that the subcells in the XFEM 
formulation are not hnite element cells. They are only used to build the XFEM quadrature 
formula. In particular, even if the subcells are close to be degenerated this does not effect the 
quality of the mesh. 

The numerical convergence results, shown in Table 1 for the case with blending cells with 
ramp correction, show a full convergence rate. As expected for bilinear hnite elements we 
observe a quadratic and linear convergence in the L 2 norm and energy norm respectively. Table 
2 shows the case without ramp correction for the blending cells. In this case, as expected, we 
observe a reduction of the convergence rate. 

6.2 Strong discontinuity 

We consider the same domains as in the case of a weak discontinuity, Di := {x G : ||x ||2 < 
0.5}, f22 := {x G : 0.5 < ||x ||2 < 1} and D := Di U D 2 . On these domains we dehne the 
following problem: 
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Figure 13: (a) Computing mesh at the coarsest level, (b) Coarsest mesh with visualized subcells; 
(c) Mesh at the second rehnement level with visualized subcells. 


DoF 

\\uh - m|| 

Conv.rate 

l|VK-«)|| 

Conv.rate 

161 

1.519e-02 

- 

7.007e-02 

- 

493 

4.026e-03 

1.92 

3.870e-02 

0.86 

1621 

9.980e-04 

2.01 

1.999e-02 

0.95 

5841 

2.533e-04 

1.98 

1.003e-02 

1.00 

21969 

6.311e-05 

2.00 

5.026e-03 

1.00 

84945 

1.575e-05 

2.00 

2.519e-03 

1.00 


Table 1: Convergence rates of 6.3 in L 2 and energy norm with conform space 14. 
Problem 6.4 (Strong discontinuity). Find the solution u 


-Am = 1 

in Qi, 

(39) 

M = 0 

on dfl, 

(40) 

Vmi ■ Hi = Ui — U2 

on F, 

(41) 

VM2 • n 2 = -Ml + M2 

on F. 

(42) 

The exact solution of problem 6.4 is the function m4 



m 4 • U 172 —t M 


(43) 

^ j 1(2 - \\x\\2) , 

X G Di, 

X G ^2- 



DoF 

\\uh - m | 

Conv.rate 

||V(«»-«)|| 

Conv.rate 

133 

1.155e-02 

- 

1.002e-02 

- 

417 

3.472e-03 

1.73 

4.222e-02 

1.25 

1469 

8.671e-04 

2.00 

2.317e-02 

0.87 

5517 

2.214e-04 

1.97 

1.176e-02 

0.98 

21293 

5.894e-05 

1.91 

6.559e-03 

0.84 

83565 

2.007e-05 

1.55 

3.986e-03 

0.72 


Table 2: Convergence rates of 6.3 in L 2 and energy norm with non-conform space 14. 
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Figure 14: Exact solution of problem 6.5. 


In Fig 14 a XFEM approximation of this problem is shown. 

Problem 6.5 (Weak formulation of strong discontinuity). Find u G Hq{Qi U 142 )? so that 

(Vn, V(p)n - (mi-M 2,95i)r - (-Ml+M2,952)r = (/, 95)0 Vyj G FTo (fli U 142). ( 44 ) 

Using the same notation as above the hnite dimensional conform subspace of i?o(fli U 142 ) 
is given by: 

Vh '■= {^h £ Hoi^i U 142) -FhlK £ Qi for K G (Tstd U %i), 

Fh\Ki ^ Qi for K G Tcut, 

GC'(f4\F)}. 

The discretized problem is then given by: 

Problem 6.6 (Discrete formulation of strong discontinuity). Find Uh G Vh, so that 

{Vuh, 'V(ph)n — {uh,! — Mh,2, Fh,i)r — {—Uh,i + Uh,2, Fh,2)r = (/, Fh)n '^Fh G 14 - ( 45 ) 


The results of the numerical convergence analysis is given in Table 3. Again, we observe 
the full convergence (quadratic in the L 2 norm and linear in the energy norm) in an unhtted 
mesh. 


DoF 

\\uh - u\\ 

Conv.rate 

l|VK-«)|| 

Conv.rate 

133 

1.663e-02 

- 

8.278e-02 

- 

417 

4.503e-03 

1.88 

4.231e-02 

0.97 

1469 

1.045e-03 

2.11 

2.155e-02 

0.97 

5517 

2.875e-04 

1.86 

1.072e-02 

1.01 

21293 

7.203e-05 

2.00 

5.364e-03 

1.00 

83565 

1.803e-05 

2.00 

2.683e-03 

1.00 


Table 3: Convergence rates of 6.6 in L 2 and energy norm. 
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7 Conclusions and possible extension of the program 

We have presented an implementation of the extended Finite Element Method (XFEM) in the 
FEM library deal.II. The implementation is mainly based on the objects hp::DoFHandler and 
FENothing. 

As part of this work, we make available a code that can be used to solve interface problems 
using the XFEM in two dimensions. The main parts of the implementation are 

• the XFEM quadrature formula; 

• the assembling routine that uses the extended part of the hnite element formulation; 

• the visualization routine for cut cells. 

We have presented two prototypical examples to numerically approximate interface problems 
with strong and weak discontinuities respectively. The numerical results show the expected 
convergence rates. In case of a weak discontinuity also the blending effect (nonconformity) and 
the resulting loss of convergence rate are shown. Furthermore, we present the implementation 
of a known remedy to restore conformity. 

We underly that a possible extension of this code is the implementation of a more efficient 
quadrature rule. The implementation shown here can be straightforwardly extended to the 
three dimensional case by dehning the necessary cell subdivisions. In practice, this leads to a 
limited efficiency of the quadrature formula. Therefore, the focus of our next work is on an 
efficient quadrature rule for the three dimensional case. 

An other extension of the program is the inclusion of extended shape functions for the 
approximation of crack propagation problems. 

A Note on how to produce the numerical results 

The source codes for two exemplary problems is included in the tar-hle xfem.tar. Extracting 
this hie two subfolders are produced, one for each problem. Note that the XFEM functions in 
the hie xfemJunctions.ee are the same for both programs. 

To get the code running, one needs to adjust the Makefile in the according subfolder. In line 36 
of the Makehle the correct path to the deal.II directory needs to be inserted. The used deal.II 
version needs to be at least version 8.0. The code is compiled with the command make and run 
with make run. 

Both programs will run performing several cycles of global mesh rehnement. 

Strong discontinuity 

The program in the subfolder strong produces the results shown in the Table 3 of the Section 
6.2. The following output is produced: 

Cycle 0: 

Number of active cells: 80 

Number of degrees of freedom: 133 
L2 error = 0.0166268 
energy error = 0.0827805 
Cycle 1: 

Number of active cells: 320 

Number of degrees of freedom: 417 
L2 error = 0.00450276 
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1280 

1469 


energy error = 0.0423088 
Cycle 2: 

Number of active cells: 

Number of degrees of freedom: 

L2 error = 0.00104475 
energy error = 0.0215455 
Cycle 3: 

Number of active cells: 5120 

Number of degrees of freedom: 5517 
L2 error = 0.000287494 
energy error = 0.0107175 
Cycle 4: 

Number of active cells: 20480 

Number of degrees of freedom: 21293 
L2 error = 7.20349e-05 
energy error = 0.0053637 
Cycle 5: 

Number of active cells: 81920 

Number of degrees of freedom: 83565 
L2 error = 1.80274e-05 
energy error = 0.00268312 
L2 Energy 

1.663e-02 - 8.278e-02 

4.503e-03 1.88 4.231e-02 0.97 
1.045e-03 2.11 2.155e-02 0.97 
2.875e-04 1.86 1.072e-02 1.01 
7.203e-05 2.00 5.364e-03 1.00 
1.803e-05 2.00 2.683e-03 1.00 

In the first part of the output information on the mesh and the errors in every rehnement cycle 
is shown. In the last part of the output a table containing the L 2 and energy errors with the 
corresponding convergence rates is reported. 

Furthermore a hie for the visualization of the solution is produced for every cycle in the vtk 
format. 

Weak discontinuity 

The program in the subfolder weak can be used to solve an interface problem with weak discon¬ 
tinuities. There is an additional parameter hie for the weak code, which contains the following 
parameters: 


set Using XFEM =true 

set blending =true 

set Number of Cycles =6 

set q_points =3 


Four diherent parameters can be adjusted. In the hrst line the XFEM is activated. Setting this 
parameter to false the program will use the standard FEM with unhtted interface. With the 
second parameter one can decide whether to use the ramp correction for the blending cells or 
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not. Remind that this parameter has no effect, if the XFEM is not used. The third and fourth 
parameter set the number of refinement cycles and the number of quadrature points for each 
quadrature formula. 

This program produces the results shown in the Table 1 of the Section 6.1. Running it with 
the parameter set as above, the following output is produced: 

Cycle 0: 

Number of active cells: 80 

Number of degrees of freedom: 161 
L2 error = 0.01519 
energy error = 0.0700709 
Cycle 1: 

Number of active cells: 320 

Number of degrees of freedom: 493 
L2 error = 0.0040257 
energy error = 0.0386981 
Cycle 2: 

Number of active cells: 1280 

Number of degrees of freedom: 1621 
L2 error = 0.00099799 
energy error = 0.0199881 
Cycle 3: 

Number of active cells: 5120 

Number of degrees of freedom: 5841 
L2 error = 0.000253274 
energy error = 0.0100287 
Cycle 4: 

Number of active cells: 20480 

Number of degrees of freedom: 21969 
L2 error = 6.31119e-05 
energy error = 0.0050257 
Cycle 5: 

Number of active cells: 81920 

Number of degrees of freedom: 84945 
L2 error = 1.57511e-05 
energy error = 0.00251584 

L2 Energy 

1.519e-02 - 7.007e-02 

4.026e-03 1.92 3.870e-02 0.86 
9.980e-04 2.01 1.999e-02 0.95 
2.533e-04 1.98 1.003e-02 1.00 
6.311e-05 2.00 5.026e-03 1.00 
1.575e-05 2.00 2.516e-03 1.00 


In the first part of the output information on the mesh and the errors in every cycle is shown. In 
the last part of the output a table containing the L 2 and energy errors with the corresponding 
convergence rates is reported. 
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Furthermore a file for the visualization of the solution is produced for every cycle in the vtk 
format. 
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